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j^ ' Abstract 

Numerical integration of ODEs by standard numerical methods reduces a contin- 
uous time problems to discrete time problems. Discrete time problems have intrinsic 

^ ■ properties that are absent in continuous time problems. As a result, numerical solu- 

tion of an ODE may demonstrate dynamical phenomena that are absent in the original 
ODE. We show that numerical integration of system with one fast rotating phase lead 
to a situation of such kind: numerical solution demonstrate phenomenon of scattering 

cn . on resonances that is absent in the original system. 
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Numerical integration of ODEs reduces a continuous time problems to discrete time prob- 
lems. For arbitrarily small time step of a numerical method such discrete time problem may 
have intrinsic properties that are absent in the original continuous time problem (see, e.g., 
[5] and references therein). In this note we describe a situation of such kind that was not 
reported before: appearance of scattering on resonances in numerical integration of systems 
with one fast rotating phase. 
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The system under consideration has the form 

(p = h g{I, f,e), y? G S^ mod 2tt. 

Here (p is an angular variable (phase), functions / and g are 27r-periodic in (p, e is a small 
positive parameter. We assume that values of function u are separated from in the domain 
U hj a. positive constant, and that functions u, /, g are real-analytic in ?7 x §^ x [0, ^o]) ^o = 
const > 0. In system ([T]) the variables J, (p are called slow and fast variables, respectively, 
the function u/e is called a frequency. System ([1]) is called a system with fast rotating phase 
[1] . Dynamics in many problems is described by systems of form ([1]) . 

Evolution of slow variables I in system ([1]) on time intervals of lengths of order 1 is approx- 
imately described by the averaged system: 

J = FiJ), FiJ) = ^J f{J,^,0)d^ (2) 



(the averaging method, see, e.g., [1]). One can use higher approximations of the averaging 
method as well \X\. We now discuss a phenomenon that appears if system ([T]) is solved 
numerically. 

Numerical integration of system ([T]) with a fixed time step n effectively introduces into 
dynamics a new "numerical" frequency Q = 27t/k. In the process of evolution of slow 
variables / the frequency uj{I)/e changes and passes through resonances with the numerical 
frequency 2tt/k,. Passage through a resonance leads to a scattering on resonance: a deviation 
of dynamics of / from the averaged dynamics described by ([2]) (see, e.g., [7], a scattering on 
resonance for the first time was discussed in [1]). This deviation depends on value of phase ip 
at the moment of passage through resonance. Consider for example, the simplest case when 
the system is integrated by the Euler method. Then the values of /, ip at the moments of 
time UK and {n + 1)k, are related as 

/„+l = /„ + /t/(/„,V?n,e), (pn+l=V^n + l^ — + l^g{In, Vn, ^) ■ (3) 

We can say that this numerical procedure integrate the time-periodic system of ODEs 
i = Kf{I,^,e) ^ 6{t-nK), (p = Ki ^g{I,Lp,e)j ^ 6{t-nK), 

n=—oo ^ ^ n=—oo 

where 6{-) is 5-function. We can use the standard identity 
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Consider Fourier series expansion for /: 

oo 
n=—oo 

Near a low order resonance ni^^^ + n2— = 0, where rii and n2 are co-prime integer numbers, 
the dynamics is approximately described by the partially averaged equations (see [7]) 

/ = F(/) + f2 Un, (/) e'^("^^+"^^*), = ^ • (4) 

q= — oo 

For such form of equations there are asymptotic formulas for amplitude of scattering on 
resonance in different situations (see [SI [7]). Assume that / is analytic in a strip |Imy9| < a. 
Then |/„| < const ■ e"*^'"'. Then under rather general assumptions the amplitude of the 
scattering is ~ a/^c"'^'"^'. For n2 = ±1 we get that the amplitude of the scattering is 
~ y/ee '^(i)'^ , Thus the considered effect is exponentially small in natural circumstances 
when K <^ e. 

We first will demonstrate the existence of scattering on resonance for iterations of maps of the 
form ([3]). To this end we will integrate numerically by the Euler method with unrealistically 
big fixed time step the following simple system of equations with one rotating phase: 

Jl = 1, /2 = COSV9, (f = . (5) 

If this system is integrated by the Euler method, then, according to formula dl]), near the 
resonance '^^^^-^ + ^2^2 = 0, the dynamics is approximately described by the system 

/i = 1, /2 = cos(v9 + n2^t), if = . (6) 

e 

Let the resonance take place at t = t*, i.e. '^ * + n2^ = 0. Denote lu = hit*), 
if^: = ip{t^). Then passage through a narrow neighbourhood of this resonance produces a 
jump in I2 given by the asymptotic formula 
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Al2 = J -— cos{^, + n2nu + sgn{ul)^) + 0{e-^) (7) 

where u'^ = w'(/i*) (c.f. f2\ for the estimate of the error term). Here u' = du/dli, we 
assume that ^'(/i*) 7^ 0. This phenomenon is called a "scattering" because the value of 
jump depends on the value of phase ip = ip^ at the moment of passage through a resonance. 
For discrete time systems these jumps were first observed in [8]. 

We take u = Ii and the following values of the time step k: e, e/2, e/5, e/10. Take e = 0.001. 
Initial conditions for each run are /i(0) = —1, /2(0) = 1, ip{0) = 0. Thus a; = Ji = t — 1, 
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(d) K = e/10 
Figure 1. Euler discretisation, I2 = cosip, u = t — 1, e = 0.001 



and at t = 1 the system passes through the actual resonance u = 0, where I2 has a jump. 
Numerical results are presented at Fig. [H One can see a jump in /2 at t = 1 due to the 
actual resonance, as well as jumps at t = 1 + 27m{e/K,), n = 1,2, due to the discretisation. 
(Plots here and in all other figures are done with MATLAB.) 



Comparison of calculation of jumps using formula 
presented in Table 1 \j. 



with numerical results in Fig. [T] is 
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A/2 (theoretical) 


A/2 (numerical) 
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47r + l 


78462.611 


0.0653 


0.0649 


8/2 


Svr + l 


315333.34308 


0.0786 


0.0782 


s/h 


207r + l 


1973427.065167 


-0.0082 


-0.0080 


£/10 


407r + l 


7895189.742154 


-0.0784 


-0.0785 



Table 1. Theoretical and numerical results for jump, Euler discretisation 



Similar scatterings on resonances occur for other maps with fast rotating phase. We demon- 
strate this for the map that appears in solving system ([5]) by 4th order Runge-Kutta method. 
We again take /i(0) = —1, /2(0) = 1, ^5(0) = 0, e = 0.001 and the same unrealistically big 
values of the time step k: e, e/2, e/b, e/lO. Results are presented in Fig. [2l One can see 
jumps in I2 at the the same moments of time as in Fig. [H 

Now let us consider scatterings on resonances for solving the following system by 4th order 
Runge-Kutta method: 



/i = l, 



COS V^ + - COS 2(^9 + ■ ■ ■ + — ^ COS lly?, = 



We again take /i(0) = -1, /2(0) = 1, (^(0) = 0, e = 0.001. Result for k = e/2 is shown in 
Fig. O We use this unrealistically big value of the time step to demonstrate existence of the 
phenomenon. 

Resonances should appear at ni-+n2— = 0, where ni and n2 are co-prime integer numbers. 
Using u = t — 1, K = 6/2, we obtain t = — — ■ 47r + 1 at resonances. 

We analyse the right hand side of I2 term by term. For the first term, cos ip, we have ni = 1, 
the resonances occur at the points t = A; ■ 47r + 1, /c = 1, 2, 3, . . . (marked as A in Fig. [3]). For 
the next term, ^ cos 2^9, we have rii = 2. The possible values of ^2 are n2 = —1, —3, —5, . . . 

(2A; - 1) ■ 27r + 1, k = 1,2,... (shown in 

Similarly, for the term ^ cos Sif, rii = 3, 

. For the term |cos47r, 

vr + 1, 37r + 1, 57r + 1, 77r + 1, . . .. 



and the time moments for resonances are t 
Fig. [3] as B, which are 27r + 1, 47r + 1, ... 

n2 = -1, -2, -4, -5, . . ., t = Ivr + 1, fyr + 1, l^vr + 1, f vr + 1 



rii = 4, n2 



-1,-3,-5,-7 



3" ' "'3 ' '3 ' '3 

and resonances are at t 



•5 



^ Value ip^, in Table [T] was calculated as follows: (/s* = (pi + e ^ J^' ui dt — (fii + 
e~^ ((<* — 1)^ — {ti — 1)^) /2, where ii is the time moment in the considered discrete grid immediately 
preceding to i*, and (pi is the value of ip for the numerical solution at t = ti. 
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(d) K = e/10 
Figure 2. Runge-Kutta discretisation, I2 = cos 99, u = t — 1, e = 0.001 
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Figure 3. /2 = cosip + | cos 2(^9 + ■ ■ ■ + ^to cos 11^9, u = t — 1, e = 0.001, k = e/2 

The corresponding points are marked in Fig. [3] as C and D, respectively. If we go on with 
this procedure, we will find every resonance for 11 terms by increasing rii from 1 to 11. 



Now let us consider the case 

1 1 

Ji = 1, I2 = cos f + - COS 2ip + - cos 3ip + 

with infinite number of terms in I2. Let us represent I2 as 



^ 



toih 



(9) 
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le^i^ + -e^'^ 



+ ■■ 
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The result of numerical simulation by 4th order Runge-Kutta method for the same initial 
data and parameters as for system ([H]) is shown in Fig. SH It is clearly seen that the 
resonances where the jump occurs are the same as for system ([H]) for the first four terms in 
I2. For later terms ^„|_i cosriiip in I2 we can find corresponding resonances, but the jumps 
become smaller and smaller due to decreasing of coefficients ^„^_i as rii increases. 

Finally, let us numerically integrate system ([9]) by 4th order Runge-Kutta method with 
the same initial conditions and parameters, as in the previous run, but with the time step 
K, = e/10 and k, = 6/20. The results are shown in Figs. |4b] and Scl One can see how 
amplitudes of jumps decay as k decays. We present a zoom of the jump at t = Svr + 1 in 
Fig. |5l This jump corresponds to the resonance 5- — — = 0. The curve in Fig. [5] looks 
"fat". The reason is that the solution is the sum of a jump curve corresponding to the term 
cosdip in the right hand side of equation (Q, and high frequency oscillations corresponding 
to, mainly, terms cos my), m = 1, 2, 3,4 , in ([9]). If we magnify Fig. |5]we would see these 
oscillations inside "fat" curve. 
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Figure 4. I2 



(c) K = e/20 
4 cos (/9 — 2 
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uj = t-l, £ = 0.001, K = e/2 
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Figure 5. I2 = , to 
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t-1, e = 0.001, K = e/20, zoom in at t = Svr + 1 



When behaviour of a system with fast rotating phase should be studied numerically, it looks 
as the most appropriate way is to use the averaging method and its higher approximations, 
like in jBj . If, however, the direct numerical integration of the original system with a standard 
constant step numerical integrator is used (e.g., to compare results with results obtained by 
the averaging) , then the deviations of the numerical solution from the exact one have form 
of jumps on resonances between the internal frequency of the system and the frequency of 
discretisation. 
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